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ABSTRACT 


This paper is concerned with the development of the theoretical 
and mathematical background pertinent to the study of steady, corotating 
solar wind structure in all three spatial dimensions (3-D) . The 
dynamical evolution of the plasma in interplanetary space (defined as 
the region beyond roughly 35 ^ where the flow is supersonic) is approxi- 
mately described by the nonlinear, single-fluid, polytropic (magneto-) 
hydrodynamic equations. We outline efficient numerical techniques for 
solving this complex system of coupled, hyperbolic partial differential 
equations. The present formulation is inviscid and non-magnetic, but 
our methods allow for the potential inclusion of both features with only 
modest modifications. We examine one simple, highly idealized, hydro- 
dynamic model stream to illustrate the fundamental processes involved 
in the 3-D dynamics of stream evolution. We find that spatial variations 
in the rotational stream interaction mechanism produce small nonradial 
flows on a global scale that lead to the transport of mass, energy, and 
momentum away from regions of relative compression and into regions of 
relative rarefaction. The magnitude of this transport is small, but 
inside 1 AU the nonradial flow can' significantly retard shock -formation 
by allowing fluid in the compressions to slip laterally, thereby 
partially relieving the stresses built up in the stream interaction. 
Comparison with simpler models demonstrates the essential nonlinear, 
multi-dimensional nature of the interplanetary dynamics. A subsequent 
paper will be devoted to the investigation of a wide range of more 
realistic model streams. 
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INTRODUCTION 

The solar corona is frequently dominated by long-lived structures 
which are organized on a global scale. Since the solar wind has its 
origin in the coronal expansion, it is natural to expect that the 
large-scale flow pattern of the solar wind reflect the three-dimensional 
(3-D) structure of the corona (.for a recent review covering the relation- 
ship between coronal structure and the solar wind, see Hundhausen, 1977) . 
Indeed, inferences drawn from a variety of spacecraft and interplanetary- 
scintillation observations provide compelling evidence that much of this 
structure survives to the earth and beyond (e.g., see Dobrowolny and 
Moreno, 1976). 

One important class of solar wind flow is the stream interaction 
for "colliding streams"), which refers to the dynamical process by which 
steady, large-scale inhomogeneities in the coronal expansion speed couple 
with solar rotation to produce significant longitudinal rearrangement 
of material in interplanetary space. The basic mechanism is illustrated 

t 

in figure 1, where we view the flow in the ecliptic plane from above the 
north pole of the sun. The fast material emanating radially from the 
.high-speed source region near the sun is subsequently aligned by solar 
rotation with the preceding slower plasma, resulting in a spiral inter- 
action structure at large heliocentric distances. Material piles up 
at the leading edge of the' stream, where the fast flow overtakes the 
slow* while material in the trailing portions of the stream is rarefied. 
Nonradial motions (aside from those Impressed upon the flow at the 
source) are driven by the lateral pressure gradients associated with 
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the compression along the interaction front. But note that the rotation 
introduces a basic anisotropy into the picture: the azimuthal gradients 

are directly driven by the stream interaction mechanism, while the 
meridional gradients (out of the plane of the paper) arise as a consequence 
of latitudinal variations in that mechanism. The latitudinal variation 
stems partly from intrinsic latitudinal variations in the corona and 
partly from the latitudinal dependence of the rotational effect. That 
is, for structures of the same spatial width, those lying nearer the 
equator suffer a greater rotational interaction than those lying nearer 
the poles. In general, therefore, we expect the meridional gradients, 
and, hence, meridional motions, to be somewhat different from their 
azimuthal counterparts. 

Theoretical interest in the corotating stream interaction has 
centered on the description of the mechanism in the neighborhood of the 
solar equatorial plane, giving rise to increasingly complex 2-D models 
whose validity depends upon the absence or negligibility .of the local 
latitudinal gradients, i.e., 3-D structure. Practical computational 
difficulties have hindered progress toward 3-D formulations of comparable 
sophistication. The 3-D linear model (S iscoe and Finley, 1-972) proves 
inadequate in the face of the large-amplitude variations in streams, 
while the kinematic description (Barouch and Burlaga, 1976) neglects 
the dynamical reaction -of the -plasma. A nonlinear 3-D approximation 
has been synthesized under .the assumption of purely radial flow 
(Suess ^t al . , 1975) , but the effects of the nonradial motions are not 
taken into account. Axisymmetric models (e.g., Siscoe and Finley, 1-969; 
Suess, 1972; Nerney and Suess, 1975) include the meriodional flow but 
ignore the essential azimuthal properties of the stream interaction. 



A nonlinear 3-D stream model could provide valuable insights on 
both observational and theoretical fronts. For example, spacecraft 
observations of the interplanetary medium are confined to a narrow zone 
about the ecliptic plane. Lest our interpretation of these data become 
ensnared in the parochial viewpoint imposed by our limited sampling 
capabilities, it is imperative that we acquire the requisite expertise 
to treat these observations in context, i.e., as nothing other than the 
ecliptic manifestation of a phenomenon that is fundamentally global in 
nature. Likewise, the analysis of interplanetary scintillation and 
radio scattering observations, which contribute useful but less detailed 
information on high-latitude and near-sun flows, would benefit from the 
availability of a more complete description of the overall solar wind 
flow structure. As a matter of theoretical interest, the strei is and 
weaknesses of established but less elaborate models could be evaluated 
and the inferences drawn from them more knowledgeably assessed. Finally, 
we could take advantage of the capabilities of a sophisticated 3-D model 
to speculate upon the properties of the solar wind flow in regions 
currently inaccessible to empirical study. 

In this paper, we develop a fully nonlinear 3— D model, of steady 
corotating streams in the solar wind. The treatment is hydrodynamic, 
polytropic, and inviscid, and describes the flow in the supersonic regime 
(beyond ~ 35 J^)-. We outline mathematical and computational techniques 
that form the basis for a family of stream models that will eventually 
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be extended to a nonlinear, 3-D MHD formulation. We consider here one 
simple example intended to facilitate the conceptualization of 3-D 
stream dynamics, contrasting the properties of our global solutions with 
those of simpler models that have appeared in the literature. A sub- 
sequent paper will be devoted to the examination of a number of more 
complex and realistic cases, which are better suited to physical 
generalization. 
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MATHEMATICAL DESCRIPTION OF THE PLOW 

Numerous published models dealing with the dynamics of corotating 
high-speed streams (e.g,, Siscoe and Finley, 1972 ; Matsuda and Sakurai, 
1972; Goldstein, 1971, Hundhausen, 1973; Nakagawa and Wellck, 1974; 

Dryer and Steinolfsson, 1976; and Goldstein and Jokipii, 1977) have been 
based upon the magnetohydrodynamic (MHD) description of the solar wind 
flow first espoused by Parker (1958). It can be shown that the kinetic 
equations reduce to the MHD formulation on the condition that we restrict 
attention to structures that are much larger than the standard plasma 
scale lengths (e.g., gyration radius, Debye length, etc.) and to regions 
some distance from the sun (r > 35 1^) where the flow is supersonic and 
the dynamics dominated by rotation and inertia (Pizzo, 1977). It is 
presumed that the solar wind can be treated as a neutral electron-proton 
gas having infinite electrical conductivity, and that heat conduction, 
viscosity, and small-scale wave phenomena may be neglected. Furthermore, 
the magnetic field may be dropped'.! rom the formulation on the grounds 
that the flow is heavily momentum-dominated" and magnetic effects therefore 
minor. This assertion is justified for the relatively broad structures 
we will consider here,' but may break down at large radial distances 
Cr > 1 AU) where the streams approach shock formation or much nearer the 
sun in the case of streams having initially sharp boundaries (Rosenbauer 
et al . , 1977; Gosling et al.; 1978; Pizzo and Rurlaga, 1977). 

The 3-D HD equations for the corotating solar wind flow are, in 
yector form in the inertial (nonrotating) frame; 
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where p Is the (proton) mass density, U the center-of-mass velocity, 

P the total isotropic gas pressure, G the gravitational constant, the 
solar mass, and y the poly tropic index'. The independent variables are 
the usual spherical-polar coordinates (r, 0, 0) , and the parameter Q is 
the equatorial angular rotation rate of the sun (differential rotation 
is neglected) . These three equations express the conservation of Mss; 
momentum, and polytropic constant respectively. 

We set y = 5/3, the' adiabatic value, to minimize the number of 
free parameters in the calculation (that is, temperature changes in the 
plasma are determined entirely hy the work performed in compressions and 
rarefactions). This choice does force the use of somewhat elevated 
temperatures near the sun if one is to match average solar wind conditions 
near the earth, leading to an overall radial acceleration' between' the 
inner boundary and 1,0 AU. However, these effects do not seriously affect 
our results, and we would rather deal with these complications than with 
the physical ambiguities introduced by a non-adiabatic y-. Finally; the 
first term in each of equations (1) arises from the relation between 
temporal and azimuthal gradients for corotating structure in the inertial 
frame, 



This relation is exact for flows that are steady in the rotating frame 
and remaiie approximately valid for slowly varying structures (i.e. , flow 



time to a given heliocentric distance small compared to the scale of 
temporal variations at the source) . 
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Equations (1) can be written in component form and manipulated 
into the pseudo-linear vector expression: 

31 = 9U = 3D £ 

3r F ’ 80 + H ’ 34> + G 
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where U is a column vector composed of the dependent variables u^, 
n, p, u^, and Uq, F and H are nonlinear matrix coefficients, and - G is a 
nonlinear five-vector driving term (the components of F, K, and G are 
listed in Appendix A) . The characteristics are readily obtained from 
the eigenvalues of H and F. We find three independent azimuthal 
characteristics 


( 2 ) 
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and three independent meridional characteristics 
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where V.- = u , - Q r sin0 

0 0 

, 2 P 

and c = y 

s • p 

The first value in each case is simply the stre aml ine of the flow in the 
rotating frame, while the next two define the envelopes of pressure signals 
from a given point. We can guarantee that the quantities under the square 
root sign remain positive and the equations hyperbolic by situating our 
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inner boundary, r , outside the sonic critical point Cat v 35 R ) » The 

r _ o > * 1 ' * >* • * » r ' • * v 0? 

integration o£ equations (?) then reduces to a relatively tractable 
initial value problem in five unknowns. 
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NUMERICAL TECHNIQUES 

Many numerical methods have been developed to cope with systems of 
coupled, hyperbolic partial differential equations such as (2) . The 
explicit Eulerian approach described by MacCormack (1971) proves 
particularly efficient for multi-dimensional problems. Our implementation 
of the technique is inviscid, which means we must forgo any but the most 
rudimentary studies of shock formation. Thus the basic formulation 
cannot be applied to the distant solar wind, where shock structures 
dominate the flow (Gosling et al. , 1976; Hundhausen and Gosling, 1976; 
Smith and Wolfe, 1976), nor to streams having sharp gradients near the 
sun and likely to steepen appreciable inside 1 AU. However, with some 
modification, it would be possible to add an explicit artificial viscosity 
and include the magnetic field in the model to permit consideration of 
these topics. 

The tedious details of our numerical method are presented in 
Appendix A; here we will discuss only the broad outlines of the technique. 
The dependent variables U = (u^., n, P, u^, Ug) are specified on a fixed 
rectangular (6, 0) mesh lying in the spherical surface, r = r Q , which 
constitutes the inner boundary. Figure 2 shows a schematic of the mesh 
superimposed on the heliocentric global coordinate system, with the 
azimuthal direction, 0, and the meridional direction (actually the 
colatitude), 0, defined as depicted. Note that the mesh covers only a 
portion of the surface, to minimize computer demands and to avoid 
numerical difficulties associated with the geometric singularities at 
the poles. Periodic boundary conditions are imposed at the azimuthal 
edge of the mesh, while the upper latitudinal boundary is a free surface. 



with the meridional dervatives there being approximated by one-sided 
differences. This treatment is permissible provided we choose the 
boundary conditions such that the stream blends smoothly into uniform 
flow near and above the upper latitudinal edge of the mesh. As long 
as this edge lies only along the outermost fringes of the stream where 
the perturbations are very small, spurious effects upon the overall 
evolution are virtually nil. The lower latitudinal boundary is handled 
the same way, unless we are dealing with a stream that is symmetric 
about the equator. In that case, the lower half of the mesh may be 
dispensed with altogether, and the symmetry conditions simplify ex- 
pression of the meridional dervatives. 

Given U ( r 0 > 6 , 0) everywhere on the mesh, we may advance the 
solution to U (r + Ar, 0, 0) by our finite-difference form of equation 
(4). The magnitude of the Ar step is, for fixed A0 and A0» given by 
the Courant condition, which requires that the values at any given point 
be causally connected to the region demarcated by the upstream 
characteristics. The Courant step-size limit for this system is given 
by 

Ar ^ max ( A r.-. -’ y-^r 

c l max j\J maxjej 

where A. and e are the characteristic eigenvalues of equation (3) . 

Should this condition be violated, spurious numerical instabilities 
arise which can quickly destroy the solution. This occurrence is easily 
avoided by updating the condition (4) every few integration steps and 
executing the computation at some prudent fraction of the specified 
limit, say, 0.‘8 Ar £ . In addition, it is necessary to check the accuracy 
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EXAMPLE - A VELOCITY-DRIVEN STREAM AT THE EQUATOR 

To illustrate the fundamental processes involved in the 3-D stream 
interaction, we consider a simple model stream in which the radial velocity- 
alone is a function of position at the inner boundary. In the real solar 
wind, we expect all flow quantities at the inner boundary to vary, leading 
to a rather complex evolution. These mixed-mode streams will be the 
subject of the next paper in this series. For this introductory study, 
however, it is imperative to work with basic, highly idealized structures 
to obtain a feel for the parameters most likely to be important in more 
realistic cases. 

For our example, we have chosen a rather simple distribution of the 
radial velocity at the inner boundary, r Q . Since the source surface is 
a sphere centered on the sun, we represent this configuration by means 
of a contour plot projected onto a globe (Figure 3) . (The projection 
is orthographic, meaning that the globe is viewed from infinity and free 
of parallax effects, though foreshortening remains.) North is at top, 
west at right — thus the contours are mapped onto the source surface 
(r Q = 35 Rq - .16 AU) just as the sun appears in the sky. For reference, 
the usual global coordinates of latitude and longitude are provided with 
a 10° spacing. 

The contours of radial velocity in Figure 3 depict our highly 
idealized stream rising from a uniform surrounding slow flow of 290 km/ sec 
to a peak velocity of 580 km/sec (marked by the "H") at the equator. 

The circular pattern and even spacing (A = 30 km/sec) of the contours 
indicates that the radial velocity variations are symmetric about the 
stream core. The overall width and amplitude of the distribution is 
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intended to mimic typical stream behavior near earth orbit (Gosling et al. , 
1972) . (A slice through this distribution at constant 0 is depicted in 

2 

Figure 4. The variations in both 0 and 0 are of the form (sin(n x)/(rr x)] , 
where |x| £ 1.) All the other variables — u^, u^, n, and f — are held 
constant on the source surface at the slow flow levels, which are adjusted 
to yield average solar wind conditions at 1 AU under the assumption of 
uniform radial adiabatic flow (Table I) . For numerical convenience, the 
boundary conditions are periodic in 0. That is, another identical stream 
is located ± 120° in azimuth from the one depicted. However, over the 
radial span of this calculation, the spacing of the streams suffices to 
guarantee that the interaction between them is negligible. 

Our specification of .the nonradial flow components at the inner 
boundary merits a brief digression here. These have been defined as 
equal to zero since it is our express intent to study the nonradial 
motions induced by the stream interaction. In fact, the neglect of 
motions which could conceivably be imposed as boundary conditions may be 
seen as only a minor handicap, since nonradial flows which are not driven 
by the dynamical interaction must decay as 1/r and thus should be 
unimportant relative to the large fluctuations that are generated by 
streams and grow with radius. 

Figure 5 depicts the resultant contours of constant radial velocity, 
density, azimuthal velocity, and meridional velocity, respectively, at 
r = 1.0 AU. (Pressure and temperature are not shown, since the initial 
conditions and the use of the poly trope law requires n, T, and p all to 
be exactly in phase at all radii.) The A below each plot specifies the 
contour spacings, with highs and lows marked by "H" and "L" . The symbol 
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T] refers to the position of the original centerline of the pattern: 
thus, the central meridian of the source surface plot (Figure 3) is 
located 50° to the right of the central meridian of the contours of 
Figure 5. This shift is due to solar rotation, as mentioned above. 

In the u_ plot (Figure 5a) we have maintained the initial contour 
spacing of 30 km/ sec to demonstrate the steepening that has occurred 
between t (Figure 2) and 1.0 AU. (The base level, of course, has risen 
because of the overall radial acceleration in the flow, as discussed 
in the mathematical section.) The high-velocity material at the core of 
the stream at r Q has overtaken the slower material to its right, resulting 
in the compression characterized by the compaction of the contours. 
Conversely, the wide spacing of the contours to the left denotes the 
rarefaction, where the fast plasma in the stream core has run away from 
the slower flow trailing the peak. Both the compression and rarefaction 
are clearly visible in the density contours (Figure 5b) . (For reference, 

_3 

observe that the density far from the stream is about 7 cm ) . The com- 
pression and the rarefaction are greatest at the equator because the 
initial amplitude of the stream maximizes there. 

Figure 5(c,d) portrays the nonradial flow components at 1.0 AU. In 
the plot, the solid contours show flow to the right (the positive or 
corotating direction) while the dashed curves indicate motion to the left 
(the negative or anti-corotating direction). Similarly, in the Uq plot 
solid contours depict positive (north-to-south) flow, while the dashed 
lines denote negative (south-to-north) motions. All of the behavior 
evident in Figure 5(c,d) is readily explained in terms of the pressure 
gradients generated by the 0-dependent stream interaction. Because the 
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stream interaction maximizes at the equator, the largest deflections 
and shears must occur there, also. However, just as the peak u^ appears 
in the region of the largest azimuthal pressure gradients, so also must 
Uq maximize along the steepest meridional gradients. Thus the peak Uq 
motions are found well away from the equator and somewhat out of phase 
with u^. In addition, u^ motions, being more directly driven, by the 
stream interaction, develop larger amplitudes and sharper longitudinal 
variations than the Uq flow. The relative amplitudes of and gradients 
in these two quantities may be adjudged by comparison of the contour 
spacings in Figure 5(c,d). (Bear in mind that the scale. A, differs 
by a factor of two in these two plots.) 

Obtaining a clear mental picture of the flow patterns implied by 
Figure 5(c,d) requires a practiced eye. Therefore, in order to elucidate 
the relation of this flow to the overall geometry of the stream, we- have 
prepared the velocity-vector plot. Figure 6. The relative magnitude and 
direction of the nonradial flow (corrected for projection effects) at a 
number of positions in the stream are indicated by the arrows. To set 
the scale, the small arrow beneath the globe has a value of 25 km/ sec, 
while the large arrow measures 400 km/sec (the minimum radial velocity 
at 1 AU, for comparison purposes) . The vectors have been superimposed 
on the density contours (equivalent to pressure contours in this velocity- 
driven case) of Figure 5b to emphasize the nonradial flow of material 
away from the compression and into the rarefaction, defining a weak 
global vortex motion. Similar large-scale vortices are also at work in 
the (r, 0) and (r, 0) planes, but their portrayal would be somewhat 
obscured by the relatively enormous magnitude of the radial velocity 


vector. 
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All the structure visible in Figures 5 and 6 — the marked asymmetry 
in the radial velocity, the large-amplitude variations in the particle 
density, and the complex tangential vortex motions— all have been spawned 
in the mechanics of the stream interaction. The radial inflow of material 
into the high-density compression region as slow plasma is overtaken by 
fast leads to a nonradial outflow on a global scale, and vice versa for 
the rarefaction region. 

The growth of this structure with radius is displayed graphically 
in Figure 7. There we view the fractional deviations 


An 

n-n 

0 

A V Au e - 

j C 

V u 0 

n 

n 

u 

u 


o 

o 

0 


of the density and nonradial flow components about n Q and u q , the local 
slow-flow density and radial : velocity, respectively. The curves in Figure 7 
refer to the peak values anywhere in the stream; An and Au^ maximize at 
the equator in this case, whereas the greatest meridional deflection 
occurs at 0 = 72° near r^, then gradually shifts toward the equator — 
lying at 0 = 77° at 1 AU — as. the stream interaction develops. It is 
evident that while large variations in the density (and pressure) are 
generated as the stream steepens, the magnitude of the nonradial motions 
remains small, even out to 1.5 AU, where the appearance of numerical 
oscillations begins to affect the accuracy of the solution seriously. 

The growth in u^ is. tempered by two. factors: first, since the interaction, 

region lies approximately along the interplanetary spiral, the orientation 
of the associated pressure gradients systematically shifts to a more 
radial alignment as heliocentric distance increases; and, second,, geometry- 
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dictates that even if the angular scale, A 0, of the gradients is kept 
constant, the spatial scale, (r sin 0) A0, must increase with r. Thus, 
in the absence of driving forces, conservation of angular momentum 
requires a decrease in the nonradial motion away from the sun, a fact 
we have exploited in setting u* and u a = 0 at r . 

Because the nonradial motions are driven by the pressure gradients 
built up in the stream interaction, Ug is sensitive to the orientation 
of the interaction front relative to the meridional plane. In our 
example, the front is exactly perpendicular to the equator at 0 = 90°; 
hence Uq = 0 there. The meridional motions that occur at higher 
latitudes are primarily edge flows, only marginally driven by the stream 
interaction. (In a subsequent paper, we will examine a case where the 
front is inclined to the equator, directly forcing material across the 
ecliptic plane.) 

The lack of meridional flow at the equator does not mean that 3-D 
effects are negligible there. On the contrary, the vortex motion in 
Figure 6 implies a transport of mass, energy, and momentum away from 
the equator in the compression region and toward the equator throughout 
the rarefaction. These motions do not, in general, balance out. We 
compute the percentage mass-flux transport at a given latitude as: 

M(r,.0) - MO:. 0) 

MC%) = M(r ~ q } x 100 (8a) 

where 

M(r, 0) = j* § d0 (r 2 pu.) 

and i denotes summation across the entire stream at fixed 0. 


(8b) 
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In Figure 8 we plot this parameter as a function of 0 at 1 AU. The 
greatest mass loss has been sustained at the equator and not at the 
latitudes where u^ is greatest. Material has t>een systematically 
transported to higher latitudes by the equatorial divergence. This 
process proceeds monotonically with radius and approaches a value of 
4.5% in the equatorial plane, at r = 1.5 ATJ, where it is gradually 
tapering off due to geometric effects. 

The mass transport parameter M (%) primarily describes the action 
of meridional flow averaged over a span of longitude and is therefore 
global in nature. A more poignant measure of these effects on a local 
scale is presented in Figure 9. The three curves refer to the maximum 
absolute value (in arbitrary units) of the radial, azimuthal, and 
meridional divergence terms in the mass conservation equation, evaluated 
at the equator as a function of radius. Clearly, beyond abou ( t 0.75 AU, 
the nonradial terms rapidly diminish in relative importance, as the 
geometric decay overpowers the continual growth in and Uq that 
accompanies the steepening of the stream. In addition, the meridional 
terms are relatively small compared to the azimuthal contributions as 
had been expected. 

These two measures of nonradial effects — the longitudinally- averaged 
quantity M (%) and the locally-sensitive comparison of divergence 
terms — together suggest that while the equatorial portion of the stream 
continues to lose mass beyond 1 AU, this divergence must be rather 
diffuse and have little local effect upon the stream evolution except 
near the sun. That is, we conclude that the observable influence of 
the 3-D geometry upon the stream development at any fixed latitude is 
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largely confined to the region inside 0.75 AU and, furthermore, that the 
evolution is more strongly affected hy the azimuthal motions than by the 
meridional motions. (The same arguments could have been applied to the 
momentum and energy divergence and transport terms; we have concentrated 
on the mass parameter because the nonradial effects are relatively 
largest in this quantity.) 

An important implication of the above findings is that this stream 
does not spread very much in latitude in the supersonic flow regime of 
interplanetary space. Further support for this claim may be obtained 
by following the motion of tracer particles embedded in the flow. 

Figure 10a depicts an array of such markers uniformly dispersed across 
the source surface, r . (The initial spacing is A0 - 2.0°, A 0 = 5.0°.) 
Differences in u^ as a function of position (see Figure 3) coupled with 
rotation lead, of course, to the dynamic stream interaction, which moves 
the particles about relative to one another . In Figure 10b we have 
tracked the markers to 1 AU, using the routine described in Appendix A. 
(This diagram has been rotated by the same angle t] as the previous contour 
plots to facilitate comparison. The lines connect the particles in their 
original sequence at r Q . ) We wish to make just two points: (1) the 

particles exhibit much greater relative displacement in 0 that in 0, as 
predicted. And (2), one may mentally overlay the u^ and n contours of 
Figure 5 to relate pictorially the compression and rarefaction to the 
local density of the convected tracers. The displacement of the markers 
in 0 is a straightforward process. However the 0-displacement of the 
markers has two components: one due to the actual motions seen in 

Figure 5c, and of a magnitude roughly comparable to the 0-displacement, 
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and' a much larger one' stemming from 1 the rotational translation of the 
original u^(0) distribution. That is-, the particles tied to the faster 
fluid elements hea r r the stream core arrive at 1 AU relatively ahead 
(to the right)' of the' markers attached to the slower preceding material. 

In the preceding' paragraphs, we have investigated the dynamical 
processes associated with a' fully nonlinear 3-D stream. It is instructive, 
therefore,- to contrast the predictions of this model with those of 
familiar, less sophisticated formulations; To this end, we direct 
attention to the evolution of the stream in the immediate vicinity of 
the solar equator; where! our 3-D solution predicts the maximum meridional 
transport effects. Figure li compares solutions for identical boundary 
conditions along the equator executed with nonlinear 3-D, 2-D* and 1-D 
models and with a 2-D linear model. Using this illustration, we will 
explore the effects attributable to the variation iii geometrical properties 
among the models, and then we will consider the role of nonlinearities 
upon the dynamics. 

Turning first to the dimensionality question, the 1-D solution 
(dashed) In Figure i 11 refers to the so-called " (r, t)" model (Hundhausen, 
1973), wherein only the radial terms in the equations (and, of course, 
rotation) .are retained. The -2-D solution (light solid) goes one step 
further and takes the azimuthal motion into account (e.g., Goldstein, 

1971) , whereas our 3-D prediction (heavy solid) includes the full 
geometry. Figure 11 displays the radial’ velocity, density, and azimuthal 
velocity variations in the ecliptic plane at 1 AU for each of the three 
geometries. They are presented in time sequence form, with time in- 
creasing to the right in the usual manner of spacecraft data display. 
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However, we have chosen to employ the more fundamental heliocentric 
angular units rather than temporal units, Fpr reference, the synodic 
angular velocity rate for an observer at 1 AU corresponds to 13.3° per 
day. Note also that time and azimuth, as defined in the spherical 
coordinate system' attached to the sun, run in opposite directions. Hence, 
this time series is a mirror image of a slice through the equator of the 
distribution in Figure 5. 

Xt is evident from Figure 11 that the 3-D geometry has had little 
influence on the evolution of the stream over that of the 2-D approxi- 
mation. The. real change comes in the transition from 2-D to 1-D (r, t) 
solutions. The azimuthal divergences have significantly broadened and 
reduced the compression at' the leading edge by allowing material to slip 
laterally away from the interaction region; the 3-D case sees a further 
augmentation of this effect. The simplified geometry of the 1-D model 
leads to a gross overestimate of the dynamical reaction to the kinematic 
steepening of the stream and to the' prediction of an exaggerated com- 
pression that is too steep and narrow. As a consequence, the 1-D radial 
velocity profile shows a stepped structure that portends shock-pair 
formation (actually occurring near 1.25 AU) . The nonradial flow in the 
multi-dimensional models, on the other hand, relieves the stresses at 
the. leading edge and thus delays the appearance of corotating shocks. 

(Due to numerical difficulties in this inviscid model, we cannot be more 
specific than to say that the multi-dimensional models -do not shock 

t 

before 1.5 AU.) We. point out, however, that once the stream reaches 
1 AU, Figure 9. suggests that the nonradial flow effects must wane, and 
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therefore the evolution at larger heliocentric distances must proceed 
in nearly 1-D fashion. 

The essential nonlinearity of the stream interaction is demonstrated 
hy comparison of the 2-D linear and nonlinear solutions found in Figure 11, 
The predictions of the linear model (dotted) were obtained by numerical 
integration of a 2-D formulation of the linear perturbation equations in 
(Siscoe and Finley, 1972) , using the computational techniques described 
in this paper. The discrepancy between the linear and nonlinear models 
are even more striking than those between the 1-D and 2-D nonlinear- 
solutions. The intrinsic inability of the linear model to steepen — due 
to the inadequate treatment of the crucial (U.y) U term in the radial 
momentum equation — precludes its application to solar wind structures 
having large variations in the radial velocity (for more details, see 
Pizzo, 1977) . On the other hand, a purely kinematic model — which may be 
viewed as the logical extreme of the 1-D nonlinear formulation — errs as 
disastrously in the opposite direction, through neglect of the dynamical 
reaction (i.e., pressure forces) of the system. (Calculations indicate 
that the kinematic model shocks near 0.8 AU for this example.) 

These comparisons imply that accurate representation of solar wind 
dynamics inside 1.0 AD demands the use of models that are both multi- 
dimensional and nonlinear. While the nonradial flow may be small 
throughout the region of interest (e.g., refer to Figure 7), their 
cumulative effect in broadening the structure at the leading edge of 
the stream may be substantial* Apparently a 2-D nonlinear model may 
suffice for the study of many aspects of stream evolution; indeed, for 
some geometries and applications, the 3-D approach may prove extravagant. 
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However the 3-D model alone provides estimates of the meridional flow 
properties and further studies will be required to distinguish when the 
effects of meridional motions may be significant. 
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SUMMARY AND CONCLUSIONS 

We have developed- a numerical three-dimensional hydrodynamic model 
to study the effects of global, corotating stream evolution in inter- 
planetary space. On a practical basis, the technique proves computationally 
efficient and is possessed of great versatility as well. With a minimum 
of effort, the model can be expanded to a 3-D' MHD formulation or, 
alternatively, reduced in scope to 2-D and even 1-D descriptions. In- 
deed, the model may be applied to any corotating flow of astrophysical 
interest, so long as the governing partial differential equations are 
hyperbolic and well-behaved. With, modest" additional labor, shocks and 
discontinuities could be accommodated through the introduction of a 
suitable artificial viscosity. 

The one simple stream model presented in this paper was selected 
to demonstrate the capabilities of our numerical techniques while also 
illustrating the basic phenomena to be encountered in corotating 3-D 
flows. A variety of graphical displays — the contour, velocity vector, 
and tracer particle plots — were utilized to elucidate the global aspects 
of the nonradial flow generated by the .3-D stream interaction. The 
primary finding is that the induced meridional motions are small, and so, 
therefore, are their effects upon the stream evolution. Thus much of 
the 3-D solar wind structure engendered in the upper reaches of the 
corona apparently persists in 1 ATI and beyond, modified in transit by 
interplanetary processes that are essentially 2-D in nature. 

Despite the hypothetical character of this example's boundary 
conditions, it is to be expected that the general sense of the results 
obtained should hold valid for a broad range of solar wind structures. 
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The flow is so heavily momentum-dominated that a significant departure 
from the traditional concepts of near-sun stream configurations— upon 
which the choice of this model was- Based— would Be required Before a 
substantial change in the dynamics - of the situation could occur. Never- 
theless, it would Be inappropriate to generalize too freely upon the 
results of this one example, and, therefore, we defer such elaborations 
to the second paper in this series. 
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APPENDIX A - DETAILS ,QF NUHERICAI frjETHQDS 
I. Differencing Scheme, and Matrices 

We wish to numerically integrate- the. 3^D HD equations in matrix 

form: 
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where the independent variables, are the usual (r, 0, 0) and the 
dependent variables are expressed as the vect.or 
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Taking equations (1) in the order radial momentum, mass, polytrope, 
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In the supersonic region, equations (A,!) are everywhere hyperholic, 
and the solution in the region G? > X ) reduces to an initial-value 
problem. That is, we are given- tj over some sp,an of 0 and 0 on the 
spherical surface r o (.see Figure -2) and are looking for a suitable 
numerical representation of equations (A.l) that will allow us to step 

I 

the solution from r to r + Ar and so on. 

'•'O'' ' ■» - 

An efficient differencing scheme suggested by MacCormack (1971) 
ajid applied by Thomas et al. (1972), is a second-order, two-step, 
predictor-corrector of the form: 

(step 1) 


«n+l _ TT n 
U. . - U: . + Ar 

30 W 


(step 2) 
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(A. 2) 

The superscript n refers to the radial mesh step, while subscript i and 
j refer to the azimuthal and meridional mesh nodes, respectively. The 
tilde denotes quantities computed from values generated in the first 
step; the vector and matrix symbols have .been eliminated for clarity, 
but it is to be understood that (A.2) represents a matrix expression. 

This differencing scheme offers one particular advantage over dozens 
of other variants that could haye been employed: storage is minimized 

since only U n (and not G n , R n , F n ) need be retained from step 1 in order 
to execute step 2. 
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Additional savings in computer time can be realized by means of the 
coordinate transform found in Goldstein and Jokipii (1976) . This 
transform is specifically tailored to- the corotating stream problem 
and secures a larger radial step size by removing most of the azimuthal 
"motion” of the stream through the fixed mesh, of the calculation. This 
apparent motion, is nothing other than the spiral geometry of the inter- 
action front. The transform effectively shifts the mesh to keep the mean 
spiral* centered* and is almost equivalent to solving the system in the 
rotating frame. 

The new independent variable in the azimuthal direction is 

§ * 0 + Qr/U r (A. 3) 

where U is the solar rotation rate and tf some average radial velocity 
at r . When the transform is applied, the radial step-size criterion (4) 
is improved since 

U1 + u/u . 

r 

Note that the exact value of is not critical; any rough average over 
the stream will serve the purpose, and, for greatest economy, need be 
updated only every few steps. One may keep track of the total absolute 
shift between the initial azimuthal coordinate system and the transformed 
system as 

A0 - - f -ii- dr 
r Q S r (r) 

The differencing, method (A. 2) can, in principle, be used to solve 
any system of equations that can be cast in the form (A.l). Thus, using 
basically the same techniques, one can tackle a variety of models — from 
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HD HD tio 3-D MHD and eyen 2— D tinie.^deLp end exit MHD— ^xi^er ely by generating 
-the appropriate, expressions, for H,* j?i and- G, The only requirement is 
that the -solutions he hyperbolic everywhere within the region of in- 
terest and that they he irell-behaved ; i.e,, hot contain discontinuities 
or gradients that afe too steep. This latter condition must remain some- 
what ill-rdefined, but; |h general, means that shocks are beyond the 

scope of the scheme described here. To Be sure, (A. 2) does contain a 

- 1 

small amount of implicit artificial viscosity and therefore may tolerate 
weak discontinuities in some cases. However for situations of greater 
interest — say, the shock-dominated structures at several AU distance 
from the sun— the introduction of explicit artificial viscosity terms 
would be required (e.g., see Ch; 6 of Roache; 1972), 

IX. Accuracy Considerations 

One of the simplest ways to test the accuracy of a numerical solution 
is to repeat the calculation using different values for* the various mesh 
parameters. This may prove convergence to a limit, though not necessarily 
to the desired one. Another method is to substitute the generated 
solutions back into the original differential equations, a ploy that is 
useful in ferreting out local errors but insensitive to systematic 
deviations that become large cumulatively. Since equations 0-) are not 
in explicit conservative form, we may resort to conservation laws to 
check our results. Easiest to implement is an integral summation of 
the conserved fluxes. If equations Cl) are cast in conservative form 
and integrated in r and 0 at fixed 0, then five relations— one -for mass, 
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total energy, and each of the momentum fluxes — can be derived in the 
form 


J d#> [flux density at (r, 0-)] + JdrJd0 
$ R $ 


change in flux density 
at 0 per unit radial 
increment 


= const 


where the subscripts 5 and R denote integration over the entire 
azimuthal extent of the calculation from radius r to r. The kernel 


o 

of the double integral contains a term arising from the exchange of 
momentum and energy among kinetic, thermal, and gravitational forms and 
a term related to the latitudinal transport contributed by the meridional 
motions. A sensitive measure of conservation on a local scale may be 
obtained through manipulation of equations (1) . Along a streamline in 
the rotating frame, we find 


-v - 

pp ' = const 

1 tn 2 . TT 2 J- TT 2 \ -L v P % -Sir U , SB 0 = const. 

2 (lJ r + B e + > + 7=1 p - — - * 

The streamlines are specified by the first characteristic in equations 
(3) and are traced with the aid of a second-order Taylor expansion about 
the mesh nodes. The expressions are tedious and convolved, and their 
integration requires an iteration loop (see Pizzo, 1977). 

These methods were used to evaluate the example presented in the 
text. For mesh resolutions A0 and A0 of 2.5° and 1.0°, respectively, 
we find 1-2% overall conservation accuracy, with 5-8% local deviations 
in steep gradient regions, rapidly deteriorating near 1.5 AU where 
the evolution proceeds vigorously. (Note that more resolution is 
usually necessary in the azimuthal direction, as those gradients are 
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more directly driven by the stream interaction. Thus A0 normally 
•dictates the- radial step criterion, A-r^ . ) The computation time for 
this example, exclusive of graphics and diagnostic routines, was approxi- 
mately two minutes on the GSFC IBM 360/91 computer. 
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FIGURE CAPTIONS 


Figure 1 Schematic of the stream interaction in the inertial frame. 

The view is from above the north pole of the sun, looking 
down on the ecliptic plane. Spatial differences in the 
nearly radial expansion (indicated by the dark vectors) 
couple with solar rotation to produce compressions (shaded) 
and rarefactions in the interplanetary medium. Secondary, 
nonradial motions are driven by pressure gradients built 
up in the stream interaction (large open arrows) . Magnetic 
field lines — which correspond to streamlines of flow in 
the rotating frame*— are drawn out into the spiral con- 
figuration as shown. 


Figure 2 Definition of coordinate, system and mesh parameters at fixed 

radius. The computational mesh (dark grid pattern) covers 
only a portion of the total spherical surface. 

Figure 3 boundary conditions for the example in the text. The 

contours of constant radial velocity are spaced at an 
interval. A, of 30 km/ sec. 

Figure 4 Longitudinal slice through the distribution of Figure 3. 

The curve is depicted as a time sequence increasing to the 
right (the usual observational convention) , while azimuth 


increases to the left. 
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Figure 5 Contours of constant radial velocity Ca) , number density 

(b), azimuthal velocity Cc) , and meridional velocity (d), 
for the stream at 1.0 Ah. The solid contours in Cc) . 
refer to westward (positive or corotating) flow, the 
dashed contours to eastward Cncgative P* ant i-corota ting) 
flow. The' solid contours in Cd) refer to southward 
(positive) flow, the dashed contours to northward 
(negative flow. The pea.lt values are. denoted by the "H ir , 
the minima by m L" , and the contour intervals by .A. The 
contours have been shifted in longitude by the angle T] 
from the original centerline to account for rotation 
Csee text) . 

Figure 6 Nonradial flow vectors at 1.0 AU, superimposed on the 

density contours of Figure 6(b), The scale is set by the 
vectors at bottom. 

Figure 7 Growth of density and nonradial velocity perturbations 

with radius. The curves refer to the peak values anywhere 
within the stream. 

Figure 8 Net percentage mass-flux transport as a function of 

colatitude, 6, at 1.0 AU. Material has been moved away 
from the equator toward higher latitudes by the meridional 
flow. 

Figure 9 Radial dependence Cat the equator) of the peak values of 

radial, azimuthal, and meridional divergence terms in 
the mass conservation equation. 
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Figure 10 Relative displacement of tracer particles Imbedded in 

flow between r and 1.0 AU. The arrow identifies the 
o 

series of particles lying along the central meridian at 
r Q . Through interpolation, one may map the interplanetary 
distortions of any pattern (e.g., magnetic topology in 
weak-field limit) frozen to the flow. 

Figure 11 Comparison of stream models in equatorial plane. Predictions 
of radial velocity, density, and azimuthal velocity at 
1.0 AU are presented for 3-D, 2-D, and 1-D nonlinear, and 2-D 
linear formulations. 
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